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Abstract 

A generalized eigenvalue algorithm for tridiagonal matrix pencils is presented. The algorithm appears as 
the time evolution equation of a nonautonomous discrete integrable system associated with a polynomial 
sequence which has some orthogonality on the support set of the zeros of the characteristic polynomial for 
a tridiagonal matrix pencil. The convergence of the algorithm is discussed by using the solution to the initial 
value problem for the corresponding discrete integrable system. 
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1. Introduction 

Applications of discrete integrable systems to numerical algorithms are important and fascinating top- 
ics. Since the end of the twentieth century, a number of relationships between classical numerical algorithms 
and integrable systems have been studied (see the review papers |H~H3"]). On this basis, new algorithms based 
on discrete integrable systems have been developed: (i) singular value algorithms for bidiagonal matrices 
based on the discrete Lotka-Volterra equation flJJ [5), (ii) Pade approximation algorithms based on the dis- 
crete relativistic Toda lattice [6 J and the discrete Schur flow [7|, (iii) eigenvalue algorithms for band matrices 
based on the discrete hungry Lotka-Volterra equation [8J and the nonautonomous discrete hungry Toda 
lattice [9J, and (iv) algorithms for computing D-optimal designs based on the nonautonomous discrete Toda 
(nd-Toda) lattice [10 j and the discrete modified KdV equation fTTTI . 

In this paper, we focus on a nonautonomous discrete integrable system called the Ru chain lfl2l , which 
is associated with the generalized eigenvalue problem for tridiagonal matrix pencils [13 1. The relationship 
between the finite Rn chain and the generalized eigenvalue problem can be understood to be an analogue 
of the connection between the finite nd-Toda lattice and the eigenvalue problem for tridiagonal matrices. In 
numerical analysis, the time evolution equation of the finite nd-Toda lattice is called the dqds (differential 
quotient difference with shifts) algorithm 1141 , which is well known as a fast and accurate iterative algorithm 
for computing eigenvalues or singular values. Therefore, it is worth to consider the application of the finite 
Rn chain to algorithms for computing generalized eigenvalues. The purpose of this paper is to construct a 
generalized eigenvalue algorithm based on the finite Rn chain and to prove the convergence of the algorithm. 
Further improvements and comparisons with traditional methods will be studied in subsequent papers. 

The nd-Toda lattice on a semi-infinite lattice or a non-periodic finite lattice has a Hankel determinant 
solution. In the background, there are monic orthogonal polynomials, which give rise to this solution; monic 
orthogonal polynomials have a determinant expression that relates to the Hankel determinant, and spectral 
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transformations for monic orthogonal polynomials give the Lax pair of the nd-Toda lattice lfl5l[T6ll . Espe- 
cially for the finite lattice case, we can easily solve the initial value problem for the nd-Toda lattice with 
the Gauss quadrature formula for monic finite orthogonal polynomials. This special property of the discrete 
integrable system allows us to analyze the behaviour of the system in detail and tells us how parameters 
should be chosen to accelerate the convergence of the dqds algorithm. We will give a review of this theory 
in Section |21 

The theory above will be extended to the Rn chain in Section [3] The three-term recurrence relation that 
monic orthogonal polynomials satisfy arises from a tridiagonal matrix. In a similar way, a tridiagonal matrix 
pencil defines a monic polynomial sequence. This polynomial sequence, called monic Ru polynomials [17|, 
possesses similar properties to monic orthogonal polynomials and their spectral transformations yield the 
monic type Rn chain. A determinant expression of the monic Rn polynomials gives a Hankel determinant 
solution and, in particular for the finite lattice case, a convergence theorem of the monic Rn chain is shown 
under an assumption. This theorem enables us to design a generalized eigenvalue algorithm. 

The dqds algorithm is a subtraction-free algorithm, i.e., the recurrence equations of the dqds algorithm 



do not contain subtraction operations except origin shifts (see Subsection 2.3 1. The subtraction-free form 
is numerically effective to avoid the loss of significant digits. In addition, there is another application of 
the subtraction-free form: ultradiscretization [18 1 or tropicalization [19]; e.g., the ultradiscretization of the 
finite nd-Toda lattice in a subtraction-free form gives a time evolution equation of the box-ball system with 
a carrier [20]. In Section |4j for the monic type Rn chain, we will present its subtraction-free form, which 
contains no subtractions except origin shifts under some conditions. It is considered that this form makes 
the computation of the proposed algorithm more accurate. At the end of the paper, numerical examples 
will be presented to confirm that the proposed algorithm computes the generalized eigenvalues of given 
tridiagonal matrix pencils fast and accurately. 

2. Monic orthogonal polynomials, nd-Toda lattice, and dqds algorithm 

First, we will review the connection between the theory of orthogonal polynomials and the nd-Toda 
lattice. 

2.1. Infinite dimensional case 

Let us consider a tridiagonal semi-infinite matrix of the form 



B <0 



wf «<*> 



w 3 
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4 6 C, ivi t] € C - {0}, 



where t e N is the discrete time, whose evolution will be introduced later. Let l n denote the identity matrix 
of order n and B^p the n-th order leading principal submatrix of B^\ We now introduce a polynomial 
sequence (x)}^ Q : 

<p^(x):=l, (p { n\x):=det{xl n -B^ ) ), n = 1,2,3,.... 

By definition, (p^ (x) is a monic polynomial of degree n. The Laplace expansion for det(xl„ + i - B^\) with 
respect to the last row yields the three-term recurrence relation 

ri t A(*) = (*-4 t) )rf ) W-4 t V?-lW' " = 0,1,2,..., (2.1) 

where we set u>q := and (p_) (x) := 0. It is well known that the three-term recurrence relation of the form 
| |2.1) gives the following classical theorem. 
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Theorem 2.1 (Favard's Theorem 112TI Chapter I, Section 4]). For the polynomials {^>n (x)}£1q satisfying the 

three-term recurrence relation ( \2.1\ and any nonzero constant h^ \ there exists a unique linear functional defined 
on the space of all polynomials such that the orthogonality relation 



C^[x m $\x)]=h\lh m , n , n = 0,1,2,..., m = 0,l,...,n, 
holds, where 

h ( n t} =h ( t) w[ t) w { 2 t) ...w$\ n = 1,2,3,..., 
and S m ,n is Kronecker delta. 

From the relation \2.2) , we readily obtain the relation 

C^<p%\x)$\x)} = h { n t} S m/n , m,n = 0,1,2 

Therefore, the polynomials {<pn\x)}™ =0 are the monic orthogonal polynomials with respect to 
Let us define the moment of order m by 

l$:=£W[x m ], m = 0,1,2 

and its Hankel determinant of order n by 
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Since the monic orthogonal polynomials with respect to C^' are uniquely determined, we then find the 
determinant expression of the polynomial 0„ (x): 
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(2.3) 



Next, we introduce the discrete time evolution into the monic orthogonal polynomials by the following 
transformation from {(pfi (x)}~ to {(p n M] 



=0" 



(x - 8^)<fi M \x) = ^l(x) + qjP$\x), n = 0,1,2, ... , 



where 



4° ==- w ,:r ... , n =0,1,2 



(2.4) 



(2.5) 



and is a parameter that is not a zero of (p„(x) for all n = 0,1,2, Suppose that {(prP ( x )}^ =0 are the 

monic orthogonal polynomials with respect to C^' and define a new linear functional £( f+1 ) by 



,(t+l)l 



[P(l)]:=£ (,) [(l-S W )P(l)] 



(2.6) 



for all polynomials P(x). Then, it is easily verified that {<pn (x)}^Lq are monic orthogonal polynomials 
with respect to CM + ^' again. Since it is shown that the monic orthogonal polynomials satisfy the three-term 
recurrence relation of the form l |2.1) , another relation 



$\x) = $ n) {x) + e ( £ ) <p^\x), n = 0,1,2 

is derived for consistency. The variable e„ satisfies the compatibility condition 

(t+ i) (t+i) (t+i) (t+ i) (t) (f) (t) 

(f+1) _ (f+1) (t+1) _ (f) (f) 
with the boundary condition 

d 



(2.7) 

(2.8a) 
(2.8b) 

(2.8c) 



The transformations \2A\ and \2.7\ are called the Christoffel transformation and the Geronimus transformation, 
respectively |22|. The discrete dynamical system \2.8\ is the semi-infinite nd-Toda lattice. 

We have seen above the derivation of the nd-Toda lattice from the theory of orthogonal polynomials. 
Using this connection, we can give an explicit solution to the semi-infinite nd-Toda lattice ( |2.8) ; the solution 
is written in terms of the moments of the monic orthogonal polynomials. The time evolution of the linear 
functional \2.6\ leads to 

Fm ^ = Fm+1 ^ ■ (2-9) 

By applying this relation to the determinant expression of the monic orthogonal polynomials ( |2.3) , the 
definition of the variable ||2.5) yields 



m _ T " T n+1 
qn ~ _(f) _(t+i) ' 

L n+1 L » 



Further, applying the orthogonality relation 1 2.2 1 to equation \2.7\ , we obtain 

r t) _ £( f )[x"4°(x)] 
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If the moments ]i„ ', the elements of the Hankel determinant T„ , are arbitrary functions satisfying the 
relation l |2.9) , then these 1 2.10) and 1 2. 11) give particular solutions to the semi -infinite nd-Toda lattice | |2.8) . 
For instance, 



(0 



r t_1 

/ x m r\(x-s {,) )cv(x)dx 



satisfies the relation < |2,9) , where CI is an interval of the real line and cv(x) is a weight function on O. If the 
integral of the right-hand side has a finite value for all m e N, then this moment gives a solution. 



2.2. Finite dimensional case 

In what follows, we shall reduce the size of the tridiagonal matrix B ^ to finite N: 
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The matrix determines a system of monic finite orthogonal polynomials. The corresponding nd-Toda lattice 
is also reduced to the case of the non-periodic finite lattice of size N: 

e i) + 4 f+1) +s (t+1) =^ ) + 4?i + « (i) ' < 2 - 12 *) 
?£? ) 4 M) = tf ) 4 t) , (2.12b) 

ej') = ejf = for all f > 0. (2.12c) 

We can solve the initial value problem for the finite nd-Toda lattice [ 2.12) through the theory of finite or- 
thogonal polynomials. 

The monic finite orthogonal polynomials {(p„ (x)}^ =0 are defined in the same way as the infinite dimen- 
sional case: 0„ (x) := det(xl n - B„ ). It should be remarked that <p^(x) is the characteristic polynomial of 
For the polynomials {<pn ( x )}n=o ant ^ any nonzero constant h!jp, there exists a unique linear functional 
such that the orthogonality relation 

C^[x m (pi i) (x)]=hi t) S mrn , n = 0,1,..., N-l, m = 0,1, (2.13a) 

and the terminating condition 

cW[x m (p$(x)] = 0, m = 0,1,2,..., (2.13b) 

hold. 

Let Xq, *!/•••/ x N-i denote the zeros of the characteristic polynomial (p]9 (x), i.e., 



<'w=n( j -4 (2.i4) 

! = 

If for simplicity we assume that these zeros are all simple, the linear functional is concretely given by 
the Gauss quadrature formula. 

Theorem 2.2 (Gauss quadrature formula 11211 Chapter I, Section 6]). Let Xo,Xi,. . .,xjv-1 be the simple zeros 
of the characteristic polynomial cpff(x). For the linear functional of the monic finite orthogonal polynomials 
{(j>n ( x )}n=0' there ex ^ some constants Cq ,Cj , . . sucn that 

N-l 

C^[P(x)]=Y l c\ t) P(x i ) (2.15) 

i=0 

holds for all polynomials P(x). Further, if 'w^' ' ', . . . , v$_-y are all real and positive, then , . . . , Cj9_. are 

also all real and positive. 

The formula [ 2.15) m eans that the monic finite orthogonal polynomials {<pn\ x )}n=o w ith the termi- 
nating condition [ 2. 13b) are orthogonal on the support set of the zeros xq, X\,..., Xpj-i of the characteristic 
polynomial (pj9 (x). 

The constants Cg , Cj , . . . , Cj^^ are calculated as 

<r = -m ' f = ' 1 N ~ 1 ' ( 2 - 16 ) 



N-l 
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where (p'^Jx) is the derivative of tp^' (x). This formula is verified as follows. Due to the Gauss quadrature 
formula 1 2.15) , the moment is given by 

, . N-l , . 

l$=tf i >[x m ] = £ cfxf, m = 0,1,2 

i=0 

This yields the relation 

N-l 

fL'+i " = E C P } (*< " x />r> 7 = 0,1 N-l, m = 0, 1,2 

z=0 
**/ 

This relation and the determinant expression of the monic orthogonal polynomials ( |2.3) lead to 

^ N-l 

= ~w n c < (*/-««■) n (*vi-xv ) 2 , y =o,i,...,n-i. 

Tj^_j i=0 0<v <i/i<N-l 
A similar calculation yields 



/ s N-l , 

T ( f ) _ rr ,(0 



\ " n c i n ( x n x v y 



i=0 0<v <V]<N-l 

Further, we have 

<Tn (*;) = I! (*j ~ x i)' 7 = 0, 1, . . . , N - 1, 

i=0 

T N-1 

These equations lead to the formula [ 2.16) . 

The spectral transformations | |2.4) and l |2.7) also work for the finite dimensional case except the Christof- 
fel transformation for n = N: 

^ +1) (x) = ^\x), (2.17) 

which is consistent with the Geronimus transformation for n = N. Equation | |2.17 l means that the char- 
acteristic polynomial <p^(x) of the tridiagonal matrix is invariant under the time evolution. In other 
words, the time evolution of the monic finite orthogonal polynomials does not change the eigenvalues of 
the tridiagonal matrix Since the time evolution of the moment is given by | |2.9[ l, we have the following 
expression for the moment: 

, N-l / f-1 \ 

(=0 \ ;=0 / 

where we define t = as the initial time. Substituting this expression of the moment [ 2.18) into the elements 
of the Hankel determinant T„ ^ and applying the Binet-Cauchy formula and the Vandermonde determinant 
formula, we obtain the expanded form of r„ : 

4° = e (nf^nW^)) n k-^ ) 2 V (2.19) 

0<r <r 1 <-<r„_ 1 <N-l \ (=0 \ ;'=0 / 0<v <v 1 <n-l / 
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Hence, we can conclude that the solution to the initial value problem for the finite nd-Toda lattice 1 2.12) is 
given by 



nit 



L n L n+1 
At) (t+1)' 
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with the expanded form of t„ < 2.19) and the expression of ci 0) < |2~T6) . 

In the rest of this subsection, we will reformulate the matrix forms of the finite nd-Toda lattice. Let 
and be bidiagonal matrices of order N: 



I 1 

(0 



,(*) 
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'N-l 
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and let <p^ (x) and ^>^/ (x) be the vectors of order N: 



(x) 



h(0 



/ \ 



Then, the three-term recurrence relation | |2.1) and the spectral transformations | |2.4) and \2.7\ are written as 



(x- S ( t ))^( t+1 )(x) = R^)^.( t )(x) + ^ ) (x), 
,p^(x) = L^^ t+1 \x). 
Note that, from < |2.14[ | and {2.21a) , 

B^<p^{x i )=x i <p {t \x i ), i = 0,1,..., N-l, 



(2.21a) 
(2.21b) 
(2.21c) 



holds. This corresponds to the fact that the zeros of the characteristic polynomials cp^ (x) are the eigenvalues 
of B'*) . Moreover, it indicates that (p^ (x ( ) is the eigenvector corresponding to the eigenvalue x ( . From [ 2~2~T) , 
we have 

x^ t+l \x) = B< f+ V t+1) « + <Pn +1) (x) 

= (L^R^+s^I^^ix) + ^ +1) (x) 

= (K( f )L« +S ( f )j N )^ t+1 )(x) + ^ +1) (x). 
This yields the matrix form of the finite nd-Toda lattice | |2.12) : 
B (t+D = L (*+D jR (*+l) + s (W)j N = R (t) L (t) + s (t)j N< 

Since l/ f ) is always regular, it is shown that the tridiagonal matrices B^' and B^ f+1 ^ are similar: 

B (M) = R (t) L (t) + s (o Jn = ( L (o y 1 ( L (o K (o + s (o Jn ) L (o = ( L (o B (o L (o. 

Therefore, the eigenvalues of B^) are conserved under the time evolution. This corresponds to the fact 1 2.17) 

that the characteristic polynomial <ftj^ (x ) is invariant under the time evolution. From the result, we can see 
that the spectral transformations < |2.21b[ | and | 2.21c[ l correspond to the LU decomposition of the tridiagonal 
matrix B^ with the shift s^. 



2.3. The dqds algorithm 

For the finite nd-Toda lattice 1 2.12) , let us introduce an auxiliary variable 



£ M) *A M) -'$v n = 0,l N-l. 



Then, equations 12.12 1 are rewritten as 
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n = 0,l,...,JV-l, 
n = 1,2,..., N-l, 

for all i > 0. 
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(2.23a) 
(2.23b) 

(2.23c) 
(2.23d) 

(2.23e) 



These recurrence equations are called the dqds algorithm. 
The spectral transformations \2A\ and \2.7) yields 

tf +1) (x)= 1 V ' "J" K , n = 0,l,...,N-l. 

X - Sw 

Hence, we obtain 



j(,t; _ T_n+i v / _ » n+i 
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with 
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By a calculation similar to the derivation of the expanded form 1 2.19) , we obtain the expanded form of o~„ 



(0 



(0 
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n ( 4 0) (*,-s( f >)n(*,-s«)i n k -% ) 2 1- 
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We henceforth assume that, for the elements of the initial tridiagonal matrix B^ , the following condi- 
tions are satisfied: Hi , ■ ■ ■ , wj^-i are ^ rea ^ w 2/ • • •' w N-i are a ^ rea ^ an< ^ positive. Then, the 
tridiagonal matrix B^ is similar to a real symmetric tridiagonal matrix. The eigenvalues xq, x\,..., x^-i of 

B( ) are thus all real and simple. In addition, the constants cP\ . . . , cl?_, are all real and positive by 
Theorem 2.2 Accordingly, the solution 1 2.20 1 and [ 2.24) with the expanded forms \2.19\ and ( 2.25) gives the 
next theorem. 

Theorem 2.3. Suppose that Mq ,m( ,...,u^\ are all real, and w^' ,^4 , ■ ■ are fl " refl ' an ^ V os ^ ve - 

Choose the parameter s^ as 

< min{xo/^i/ • ■ forallt>0. (2.26) 

Then, the variables q„, e\P and d\P of the dqds algorithm 1 2.23) are positive for all n and t > 0. 



This theorem guarantees that, under the assumption, the dqds algorithm does not contain subtraction 
operations except the parameter terms -(s( f+1 ) - s^)) in equations [ 2.23a) and < 2.23b) . Namely, equations 
( 2T23J are the subtraction-free form of the finite nd-Toda lattice (2.12 1. It is known that this form improves 
the accuracy of the numerical computation. 

The asymptotic analysis of the dqds algorithm proves convergence of the algorithm and provides a 
method for accelerating the convergence. The solution derived in Subsection 2.2 is a key tool for the analysis. 
Arrange the eigenvalues Xq, x\,. . of the initial tridiagonal matrix in descending order: Xq > x\ > 

■■■> 11 trie parameter s^ f ^ chosen as 1 2.26) , namely s^ f ^ < x^-i, then the inequality xq - > x\ - s^' > 
■■■ > xn-i - s^' > holds for all t > 0. Under this assumption, by the solution | |2.20[ | with the expanded form 
[ 2.19) , we obtain the asymptotic behaviour for t » 0: 
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This shows that q„ and converge to x n -s^) and as t -* +00, respectively. Hence, it is shown that the 



dqds algorithm ( 2.23) with appropriate parameters s^ f ^ computes the eigenvalues of a given real symmetric 



=(') 



tridiagonal matrix. It is clear that the convergence speed depends on — " h (t) . Therefore, we should choose 

the parameter < x^-i as close as possible to the minimum eigenvalue for fast computation. The 
acceleration parameter s'*' is called the origin shift. 



3. Rn polynomials, Rn chain, and generalized eigenvalue algorithm 

We shall extend the discussion in Section [2] for tridiagonal matrix pencils and its associated nonau- 
tonomous discrete integrable system. 



3.1. Infinite dimensional case 

Let us consider two tridiagonal semi-infinite matrices in the following forms: 
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Let A„ and B„ denote the n-th order leading principal subma trices of A^) and B^), respectively. We now 
define a polynomial sequence {cpn ( x )}^=o by 

<pj f) (x):=l, cpi t) (x):=det(xBi t) -A { n t} ), n = 1,2,3,.... 

The polynomial cpn i x ) is a monic polynomial of degree n. In the same manner as in the case of monic 
orthogonal polynomials in Section|5J we obtain the three-term recurrence relation 



fn+i (x) = (u n t) x- v^p ) <p%> (x)~ ufp (x - K t+n _ x ){x-A n ) (p^\ (x), n = 0, 1, 2, . 



(3.1) 



where we set w^ := and cp^(x) := 0. We will assume in what follows that all the parameters K t+ / C and A/ c , 
k = 0, 1,2, . . ., are not zeros of the polynomial fn'(x) for all n e N. The polynomials {(pn\x)}^ =0 are called 
the Rjj polynomials with respect to CP', introduced by Ismail and Masson flTl . 
We introduce the notations 



K^\x):=Y\{x-K t+j ), K k {x):=K^{x) = Y\{x- Kj ), L ; (x) : =n(*-Ay), 



k-l 



,(0) 
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and V(C^ ) a linear space spanned by the rational functions 



,M = 0,1,2,. ..;m = 0,l,...,fc + Z. 



The following Favard type theorem is proved. 

Theorem 3.1 (Favard type theorem for the Rn polynomials [17]). For the Rjj polynomials {cp n ^ (x)}£L an ^ an V 

nonzero constants and h^', which satisfy * fej , there exists a unique linear functional defined on V{C^) 
such that the orthogonality relation 



£ (0 



x m cpi t) (x) 
K^\x)L n {x) 



/in <W, 71 = 0,1,2,..., m = 0,l,2,...,n, 



holds, where h„ J , n = 2, 3, . . . , are some nonzero constants. 

In the rest of this paper, we consider the monk Rn polynomials, i.e., the case where w„ = 1 + it>„ holds 
for all n = 0, 1,2, ... . For general tridiagonal semi-infinite matrices of the form if det(B^ ) + holds 

for all n = 1,2,3, . . . , then { 1 

is valid for such matrices. 

The moment of the Rn linear functional is introduced by 



are the monic Rn polynomials. Therefore, the following argument 



u k,l,t ._ At) 
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(3.2) 



and its Hankel determinant by 
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Then, the determinant expression of the monic Rn polynomials {<p B (*)},^lo * s presented 
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(3.3) 



The discrete time evolution for the monic Rn polynomials is introduced by an analogue of the spectral 
transformations for monic orthogonal polynomials: 



(x - s (t) )(l + qi t} )cp n t+1) (x) = f^(x) + qi' } (x - K t+n )<PnHx) r 
(1 + e^)qfp{x) = (p ( „ t+1 \x) + e ( „\x - A n )^_ + 1 1) (x) 
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(3.4a) 
(3.4b) 



for n = 0, 1,2, . . . , where 



Ji) (M)\ 

4° :=-(sW- Kt+n y n £\ \ n = 0,1,2 

and e„ is the variable determined by the compatibility condition: 



1 + w 



(t+i) 



-4 f+1) -4 f+1)1 ^C! + (i + 4 f+1) )(i + 4 f+1) ) 



1 (f+i) 



1" —^Ty " e „+i + (1 + ^0(1 + Ci)' 



(0 



,(f+l) 



1 (f+1) 



I (t+l) 

1 + r n -i 



-I 

-Kt+nq n ~ (J) ~ A n+l e n+ i + s ( 1 + t ?n H 1+e „ + lJ' 



(0 



1 + e 



(t+i) 



(t+l) _ ( t+ lWt+l)l+jj _ „(tWf) 



1 + e 



(0 
n+i 



fn-l e " 1 (t+l) e " (0 



with the boundary condition 
e { J } = for all f > 0. 



(3.5) 



(3.6a) 



(3.6b) 
(3.6c) 

(3.6d) 



This is the semi-infinite monic type Rji chain. Note that, since [ 3.6a) and | 3.6c| are identical, there are the two 



independent equations that determine the time evolution of the two variables <j„ and e„ ' . It is readily veri- 
fied that if {(Pn\x)}™ =0 are the monic Rjj polynomials with respect to then polynomials {cp n ^ { x )}^Lq 
defined by the spectral transformation | |3.4a) are again the monic Rn polynomials, where the corresponding 
Rn linear functional is defined by 



£( t+1 )[R(x)]:=£« 



x - s 



(0 



X-Kt 



-R(x) 



(3.7) 



for allR(x) eV(C^). 

Let us derive a solution to the monic type Rn chain. By the definition of the moment < |3.2| and the time 
evolution of the linear functional \3.7), we obtain the relations 

(3.8a) 
(3.8b) 



,,k,l,t ,,fc+l,U v „k+V,t ,,k,l+l,t j ,XM,t 
Fm ~ F m +1 ~ K t+k} l m ~ Fm+l ~ A l+\Vm 



k,l,M _ u k+V,t_Jt) fc+Ut 
rm ~ rm+i s rm 



The relation 1 3.8b) , the determinant expression of the monic Rn polynomials \33\ , and the definition of the 
variable |3.5 1 lead to 



n,n,t n,n+l,t+l 

^M^-wr 1 n n+1 



n-l,n,t+l n+l,n+l,f ' 
T " T n+1 



(3.9) 



Next, the relation < 3.8a) and the spectral transformation < 3.4a| yield 



+ <?« - <K f+ „ s J (f+1) -I s K f+»J «+l,n+l,i «-i,«,t+r 



0« {Kt+n) 



n+1 
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Further, the Jacobi identity for determinants [23 Section 2.6] proves the bilinear equation 



k-l,l-l,t k,l,t _ k-\,l,t k,l-l,t _ k-l,l-l,tjc,l,t 
'■n '■n '■n tj n-1 



By using this bilinear equation and the three-term recurrence relation l j3.ll , we obtain 



w 



(0 



£(0 








K^(x)L n+1 (x) 


4° (*)!«(*). 


n- 


£( f ) 




X»-VW(.T) 


t" 


J$\x)L n {x) 







T n+1 ' T n+\ 



Hence, from equation 1 3.6c) and these formulae, we find a solution 



(t) _ T» n 1+ Vl _ / B (t) „ \ ^n-l ^+1 



^n-\,n-\.t+\n+\,n,t 



(0 1 (0 
<7„_i 1 + <7it 



(3.10) 



If the moments are arbitrary functions satisfying the relations < |3.8| , e.g., 



# = r 

n Jn 



K t+n (x)L„(x) 



-co{x) dx, 



then ^9) and pip) give a solution to the monic type Rn chain ( |3.6) expressed by the Hankel determinant 

i n 

The reason why the Hankel determinant appears can be explained from the point of view of the dis- 
crete two-dimensional Toda hierarchy [24 j. Note that there is another determinant expression of the Rn 
polynomials and a solution to the Rn chain: the Casorati-type determinant solution I25ll26l . 

3.2. Finite dimensional case 

In this subsection, we will derive the solution to the initial value problem and the convergence theorem 
for the monic type finite Rn chain. 

Let us start with a pair of tridiagonal matrices of order N: 



A (t) 



I v {t) 
i (0 



Kt 

(0 
1 



Kf+1 



, (f) 
A2W2 



\ I 1 



Kt+N-2 
V N-1 I 



w\' 1 + w \ ' 



III 



1 

w 
(0 



\ 



(0 1 (0 



(3.11) 



The corresponding monic type finite Rn chain is 



Kt+ n +iq n + A n e n rmT _s + n i + e n ) 

l + (? i 

' n-l 



a (t) 



1 + e 



(0 



1 + e 



-^- + A„ + ie n+1 -s w (l+^„ )(l + e, !+1 ), 



(t+l) (t+l) 1 + _ (t) (0 1+g n+l 

"Jn-l e " 1 (t+l) fw e " 1 (f) ' 
1 + <?n-l 1+e n 

e W = e g) = for all t > 0. 



(3.12a) 

(3.12b) 
(3.12c) 
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To derive the solution to the initial value problem for the monic type finite Rn chain ( |3.12) , we consider 

the monic finite Rn polynomials {fn\ x )}^Lo defined by fn(x) ■= det(xB^ - An). We should remark 

that tp^(x) is the characteristic polynomial of the tridiagonal matrix pencil (A^\B^); the zeros of the 

polynomial <pw (x) are the generalized eigenvalues of the matrix pencil (A^,B^ ), i.e., the solutions of the 
equation 



A«<J> = xB (t) 0, xeC, OeC N -{0}. 
Let 2?(£W) be a linear space spanned by the rational functions 



, m = 0, 1, 2, For the monic 



finite Rn polynomials (x)}^ =0 and any nonzero constant there exists a unique linear functional 
defined on T>(C^ ) such that the orthogonality relation 



£ (0 



x m cpi t) (x) 
4° (*)**(*) 



K> 5 mA , n = 0,l,...,N-l, m = 0,l,...,n, 



(3.13a) 



and the terminating condition 



£ (0 



L K^\x)U{x) 



0, 



fc,/ = 0, 1,..., N, m = 0,1,2,..., 



(3.13b) 



hold, where the constants /Iq , ftj , . . . , are given by solving the following linear equation 



/ -1 

(0 



-(l + wf) 1 



(0 



i -(i+wJ5-iV 



"q 





\ o / 



i.e., 



, (f) rrfO/'l (0 (0 (0 (0 (0 (0 \ 

, (t) (t) (f) (f) (t) (0 (t) \ 



Note that, for the infinite dime nsional case, there are two degrees of freedom: the choice of the two constants 
/Iq and hff (see Theorem 3. 11. For the finite dimensional case, however, there is only one degree of freedom: 
the choice of the constant The cause of this is the terminating condition f3.13b[ |. 

To derive a realization of we give a quadrature formula for the Rn linear functional. Suppose that 
all the zeros xq, x\,..., x^-i of the characteristic polynomial <p(/ (x) are simple. 

Theorem 3.2 (The quadrature formula for the Rn linear functional). Let xq, X\, . . . , xn_i be the simple zeros of 
the characteristic polynomial q)j9(x). For the linear functional of the monic finite Ru polynomials {cpn\ x )}n=o> 
there exist some constants cj, , c^\ c^_ x such that 



N-l 

£<*>[*(*)] = Zcf>R(xi) 

i=0 

holds for all R(x) e £>(£«). 



(3.14) 
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Proof. This proof is an analogue of the proof to the Gauss quadrature formula (Theorem 2.2 1. For the given 
rational function R(x), consider the following interpolation rational function 

N-l 

A(*):= E 4 ° (*)*(*«), 

i=0 



where 



^(x):= 



(x - x,)^ (x)L N (x)</^ } (x ; ) 



f = 0,1,...,N-1. 



It is readily shown that 

ef\xj)=S ijt i,j = 0,1,..., N-l, 
holds. Let 

Q(x) :=R(*)-A(x). 

Then, the numerator of Q(x) is a polynomial that has zeros at xq, x\,..., x^-\- Since R(x) e V{£^), there 
exists a polynomial P(x) such that 



P(x)<P N } (x) 

K^(x)L N (x) 



By the terminating condition < |3.13b) , we obtain 

£^[R(x)] = £^[A(x)] + £^[Q(x)] 

N-l 

= ££(%f ) (x)]*(x,) + £ (t) 



!=0 



P(*)pff(*) 
L 4 } (x)L N (x) 



= ££«[4 f) (x)]R(x t ). 

! = 

Set c[ f) := £<*) [^ f) (x)], z = 0, 1, . . . , N - 1, then the proof is completed. 



□ 



Zhedanov [13 j derived a formula to calculate the constants Cq\c^\ c^_ v He used the second kind 
polynomials to derive it. Here, we give a direct calculation to check his result. From the quadrature formula 
| |3.14[ , the moment is written as 



u m - £ 



(0 



K^\x)L } {x) 



c (t) x m 
i i 



N-l 

E 

j=o KY> (xfiLiixi) 



In the same manner as in Subsection 2.2 we thus obtain the following formulae for j = 0, 1, . . . , N - V. 



9 { N-i( x j) 



1 N-l cf\x r Xi) 

N-l,N-l,t 1 1 



l ;V 1 T=0 X^ 1 (X/)L N _ 1 (X / ) 0<i/ <V<N-l 

i*i vo*j' v i*i 



n ( Xv i xv o) ' 



(t) n_i 

? N (*;) = Yl( x j- X i)' 

i=0 
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and 



N-l,N-l,f 



N-l 

n 



i=0 ^N-l( X ') L N-l( x i) 0<v <V<N-l 



( Xv i Xy o) ' 



"N-l 



N-l,N-l,t 
T _N 

N-l,N-l,t ' 
r N-l 



Hence, we find the formula 



(f) _ ^N-l^N-l( X ') L M-l( X ') 



z = 0,l,...,N-l. 



For the finite dimensional case, in the same manner as for the monic finite orthogonal polynomials (see 
Subsection 2.2 \, the characteristic polynomial is invariant under the time evolution: 

From the results in Subsection |3.1[ we can thus see that the solution to the initial value problem for the 
monic type finite Rn chain is given by 



—n,n,t n,n+\,t+l 

^(sW-xfc*)- 1 " T " +1 



n-l,n,f+l ji+l,m-l,f ' 



n-\,n-\,l+\ n+\,n,t 
,(0 _ (At) „ \ Vl T n+1 



(s W -K f+n )- 



n,n—\,t n,n,t+l 



where, because the moment is concretely given by 

i=0 Kt+k( x i)Ll( x i) 
the expanded form of the Hankel determinant is 



jy,t 



E 

0<r <r 1 <-<r„_ 1 <N-l 



'nf§<^ n ^ 



The solution derived above yields the following theorem. 

Theorem 3.3 (Convergence theorem for the monic type finite Rn chain). Suppose that all the generalized eigen- 
values Xq, %\, . . . , of the initial tridiagonal matrix pencil ( A^°), B^°) ) are real, simple and arranged in descending 

order as xq > x\ > ■ ■ ■ > xn-1- Choose the parameters s^) and Kf+N-i as x N-i > s ' f ^ an d *n-i » K t+N-if° r t >0, 
respectively. Then, we have the asymptotics of the variables for t » 0: 



(0 



s(0 



K t+n 



+ O I max • 



nj =0 (*« - s(/) ) WjTHxn-i - Kj) nj =0 (*« + i - s0) ) nps~\*n - */) 



n$(x n -i - s(f)) n^-Hxn - Kj ) ' n$(*n - so - )) u^~\x n+1 - Kj ) j ) ' 



2 " = \n; =0 (^-i- s o)) nfZixn-Kj) i 



Hence, the variables and e\P converge to ^ | - and as t -*■ +oo, respectively. 
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This theorem implies that, from l |3.6) , the elements v„ and of the tridiagonal matrices A^ and 

B (t) 

converge to x n and as t -*■ +00, respectively. Further, we can see that the parameters s^ f ^ and Kf +n 
determine the convergence speed; the parameter s^' works as the origin shift, which is the same as for the 
dqds algorithm (see the end of Subsection |Z3|. 

Next, we discuss the matrix form of the monic type finite Rn chain. Introduce the rational functions 
defined by the following three-term recurrence relation: 

*i?(*):=0, *J f) (ac):=l, 

(x - K t+n )<$> { 2 r (x) := - ((1 + w { n f) )x - pW)*W(x) - wi*\x - AOfJ^i), n = 0,1,.. .,N- 1. (3.15) 
By comparing to the three-term recurrence relation \3.1\ , the relation 

<PnH*) 



is verified. Let 



* (f )(x) := 



n = 0,l,...,N, 



K )(x) l 




f \ 














Then, the three-term recurrence relation 1 3.15) is rewritten as 
Further, let LjP, Lg , and 2?^ be bidiagonal matrices: 



(0 



-X t e\ 



(0 



1 -(t) 

-A]V-lC N _! 



/ 1 



(0 



K t+N-l) 



(0 



R (0 



(0 



(3.16a) 



.(0 

-N- 



1 V 



and Dq f ^, Dg and Dg^ be diagonal matrices: 

D? := diag (l + q^A + q^ 1 + &), 

D< f) := diag (l,l + 4° 1+4^), £>e° ^diagfl+e^ l + eg^l). 

Then, the spectral transformations | |3.4[ | are written in terms of the rational functions ( x )}^ = o as 

(x - s (t ))Dj O (t+1 )(x) = (x - Kt) (i* (f) * (f) (x) - <&j^(x)) , (3.16b) 
D< f) *<*>(*) = (x - Kt)- 1 (*4° - <5( t+1 )(x). (3.16c) 
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Equations l |3.16[ l yield 



x ( B (W)0( t+1 )(x) + < +1) (x))-A( f+1 )a»( t+1 )(x)- Kt+N 0^ +1) (x) 
x((-D^K(0 (D (0 ) -i L (0 +D (0D(0) <& ( t+ i) (x) + o a + i) (x) ) 

-(-£><Pr«\d<P)- 1 lV + ,Vd^ 



Hence, the compatibility condition for [ 3.16) , i.e. the matrix form of the monic type finite Rn chain, is given 
by 

A (M) = _ D (W) L f 1)( D (W))-1 R (W) +S (M) D (M) D (M) 
B (t + i) = _ D 4 t+1 >4 f+1 > (D4 t+1 > + Dj t+1) D( f+1) 

This leads to 

A (*+D = £)( f ) i? ( f ) (D J) D (0 ) -i A ( t ) (i? (t) ) -i D J) / 
B (*+D = DWR(') (D (0 D (0 ) -i B (t)( R (0)-i D (0 / 

and 

xg(f+1) _ A(t+1) = £(0 R ( t ) (D (0 D «)-i ( xB (0 _ A (0) (j r(0)-i D W. 

The last equation implies that the generalized eigenvalues of the tridiagonal matrix pencil (A^\B^) are 
conserved under the time evolution. 

4. Generalized eigenvalue algorithm 

4.1. Subtraction-free form of the monic type Ru chain 

In Section [5J we have presented the convergence theorem for the monic type finite Rn chain (Theo- 



rem 



3.3 1. This theorem allows us to design a generalized eigenvalue algorithm for tridiagonal matrix pencils. 



Consider a pair of tridiagonal matrices of order N as input: 



A : 



/ a 0,0 a 0,l 
01,0 0-1,1 
«2,1 

\ 



«1,2 



fl N-l,N- 



«N-2,N-1 
■2 a N-l,N-lf 



( b 0,0 



oo.i 
Ki 
h,i 



h,2 



°N-2,N-1 
V-l,N-2 b N-l,N-l! 



(4.1) 



Suppose that all the subdiagonal elements &o,i> ^1,2/ ■ ■ ■ > ^N-2,N-1 an d ^1,0' ^2,0/ • ■ ■ > ^n-i,n-2 or tne matrix B 
are nonzero, and all the leading principal minors of the matrix B are nonzero. Then, the transformation 



A (0) := ViUAirVa, B (0) : = V x UBU~ l V 2 
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gives the initial matrix pencil of the form 1 3.11) for the monic type finite Rn chain, where 
U := diag(l, b 0/ i, b 0) ib 1>2 , b 0)X b li2 ■ ■ ■ b N _ 2 ,N-i), 

Vi := diag ( (det B x ) ~\ (det B 2 )~\ . . . , (det B N )" 1 ) , V 2 := diag( 1, det B x , det B 2 , . . . , det B N _i ), 

and B n is the n-th order leading principal submatrix of the matrix B. Namely, the elements of A^ and B^ 
are computed by 

(0) ._ detB„ (o) detB„_i fl„, n+ i fl„,„_i , , 

v n — a n,H j . D / ^>n — Vn-l,nVn,n-l j . D ' Kjl ,_ 7 / A n — 7 ■ (4./) 

aetc n+1 detis„ + i f«,n+l °n,n-l 

Note that, if n is large, an overflow may occur when one computes det B n directly. The values g^fr^ and 
g"~j should be computed by the LU decomposition. Next, by the relation l |3.6) , "decompose" the matrix 
pencil (A( \ b' ') to the variables of the monic type finite Rn chain: 

4° } :=0, e^ 0) :=0, (4.3a) 
(0) 1 (0) 

4 0) :=^ 4 0) -4 0) ^g, n-1,2 N-l, (4.3b) 

rm v {0) - s(°) CI + w (0) ) - (s(°) - A„)e (0) 
4 0) :=^ { " ' [ n = 0,l N-l. (4.3c) 

Notice that the initial matrix pencil (A^°\ B^ ) does not fix the values of the parameters s^ ^ and Kjyj-i- We 
must choose the parameters s^ ) and kj^-i appropriately. We will discuss how to choose the parameters in 
the end of this subsection. After that, compute the time evolution of the monic type finite Rn chain by using 
1 3.12| iteratively; i.e., for each t > 0, compute 

4 t+1) :=0, ^-O, (4.4a) 
(t) 1 1 (0 

(M). (Q qV 1 + 1 + e „ + i 12 Ar _ 1 (44b) 

♦ (s'^-A^W -(s^ 1 )-sW)(i + ,(*))(i + e (0 1 )J / 

n = 0,1,..., N-l. (4.4c) 

Here, we also have to choose the parameters s^ f+1 ^ and Kf+jv for computing the above recurrence equations. 
From the results in Subsection |3.2| we can see that if the absolute values of all the subdiagonal elements 

A n Wn an d v>n of the matrix pencil (A^\B^ ) become sufficiently small at a time t, then the values (s^ - 
Kt+ n )cjn + s ^ §i ye th e generalized eigenvalues of the initial tridiagonal matrix pencil (A, B). In general, 
however, equation 1 4.4c) requires subtraction operations, which may degrade the accuracy by the loss of 
significant digits. A subtraction-free form of the monic type finite Rn chain may resolve the problem. 
Let us introduce an auxiliary variable 



V s K t+n+l) c ln ~ V s A n+l) e n+ i 



,(t+i)_ : ^ ■ 111L j 1±L/ n = 0,1,..., N-l. 



1 (0 
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This is an analogue of the auxiliary variable 1 2.22) introduced in the dqds algorithm. Then, the subtraction- 
free form is derived as 



l 



In 



(0 



( S ( f+1 )-s«), 



,(f+i) 



(0 

,(t+i) % 



n-l (f+i) 
ln-1 



(s< w >-s«)(l + tf>), 



n = l,2 N-l, 



0, 



S(*+D-K t+B+ 1 

,(0 

( i ! i 

e 



n = 0,l N-l, 



ii i i (f+i) 1 CO 

„(*+!)._ .(0 1 + 1 + g n + l 



n = l,2,...,N-l, 



0. 



(4.5a) 
(4.5b) 

(4.5c) 



From the spectral transformations | |3.4) , we have 

At) f^^h rrxfi +w 

where 



n+l,n+l,t-l ^_n-l,n,t+l ' 
r n+l 



s P'i+i \0<i,j<n-l- 



In addition, we already have the expression of q„ \3.5\ . Hence, we obtain a sufficient condition for com- 
puting the recurrence equation ( |4.5) without subtraction operations except the shift terms in 1 4.5a I: for all 
n and f, 



w 



(0) 



>0, 



(4.6a) 
(4.6b) 
(4.6c) 
(4.6d) 

By l |4,2) , if the input tridiagonal matrix B is a real symmetric positive (or negative) definite matrix, then the 
condition (4.6a I is satisfied. Further, assume that the generalized eigenvalues Xq,X\,. . .,1^-1 of the input 
tridiagonal matrix pencil (A, B) are all real and simple, the matrix A is a real matrix and the conditions 



(-l)>«(s«) = det(4 t) -s (t) 4°)>0, 
(-l)>«(s(^)) = det(A( t )-s( f+1 )B«)>0, 
S ( f )>A„. 



S (t) > K t +n, 



(0 



> and Kf+ n -i = A n are satisfied for n = 1, 2, . . . , N - 1 at some time t. Then, it is shown that if the 
parameter s(0 is chosen as s W < minjxo/ X\,..., *N-l}/ the condition | 4.6b[ is satisfied. The condition ( |4.6c) 
is also satisfied with s^ f+1 ^ < min{xo,xi, . . From Theorem 3.3 if s'O is chosen as close as possible 

to min{xo,Xi, ■ ■ ■ , x N-l} under the conditions l |4.6t , the convergence speed is accelerated. 

By summarizing this subsection, Algorithm [T] is proposed as a new generalized eigenvalue algorithm 
for tridiagonal matrix pencils based on the monic type finite Rn chain. 

4.2. Numerical examples 

We shall give numerical examples. To construct test problems with known generalized eigenvalues, let 
us consider the monic finite orthogonal polynomials {pn(x)}^ =0 defined by 



'n+l 



I N-l\ n(N-n) 

(x) := Ix ^—JPn(x) Pn-l(x), tl = 0,1, ... ,N -1, 
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Algorithm 1 The proposed generalized eigenvalue algorithm based on the monic type finite Rn chain 



function GEVRII(A,B) 



Compute {vi 0) }l- \ {w^J, { Kn }^, and {A„}^ by gT2 
Set the parameters and k^_i appropriately 



> A and B are tridiagonal matrices of the form ( |4.1) 

and the condition (14. 6b 



Compute {qPy^o 1 and {4 0) }£L by g3} 
f :=0 
repeat 

Set the parameters s' t+1 ^ and K f+ ^ appropriately 

Compute and {e< t+1 > by gS) 

f := f + 1 

for n = 1,2,..., AT- 1 do 

,(0 ._ „(*) „(*) 



'n— 1 



end for 

until the absolute values of w„ and k n w\p 

return {(s^-Kt+^^+sW}^ 1 
end function 



D> See Theorem 



3.3 



D> See Theorem 



3.3 



and the condition i 



are sufficiently small for all n = 1, 2, . . . , N - 1 



with p_i(x) := and po( x ) := 1- The polynomials {p»(*)}^lo are the monic Krawtchouk polynomials with a 
special parameter and it is well known that the Krawtchouk polynomials are orthogonal on x = 0, 1, . . . , N - 1 
with respect to the binomial distribution [27J. This means that the tridiagonal matrix of order N 

/(N-l)/2 1 
(N-l)/4 (N-l)/2 1 

2(N-2)/4 (N-l)/2 1 
3(N-3)/4 ■•. 



N •= 



NxN 



V (N-l)/4 (N-l)/2/ 

has the eigenvalues 0, 1, . . . , N - 1. The symmetric tridiagonal matrix 

/ (N-l)/2 ^(N-l)/A 
\/(N-l)/4 (N-l)/2 \/ 2 (N-2)/4 

v / 2(N-2)/4 (N-l)/2 y/3(N-3)/4 
V3(N-3)/4 



N : = 



V(N-l)/4 
\/(N-l)/4 (N-l)/2) 



eR 



NxN 



is similar to Kjy. Hence, it is readily shown that the tridiagonal matrix pencil (JCjy + 21^, Kj^ + 1^) has the 
generalized eigenvalues (n + l)/n,n = 1,2,. ..,N. 

The following experiments were run on a Linux PC with kernel 3.7.4 and gcc 4.7.2 on Intel Core i5 
760 2.80 GHz CPU and 4 GB memory. All the computations were performed in double precision and the 

stopping criterion (line 13 in AlgorithmjlJ was \w\P \ < 10~ 20 and \X n v$ \ < 10~ 20 for all n = 1,2, . . . ,N — 1. 
Example 1. The first example is the case of N = 5: 

\ 
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12 1 






1 2 








2 






n/3/2 


2 


\ 




1 
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2/ 
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Figure 1: The behaviour of the variables of the monic type finite Rn chain for the input tridiagonal matrix pencil (K5 + 21$, K$ + I5) 
with the parameters = 1.01 for all f > and k„ = 1 for all n > 4. 



The generalized eigenvalues of the matrix pencil (K5 + 2/5, K5 + I5) are 2, 3/2, 4/3, 5/4, and 6/5. By this 
example, we will observe the behaviour of the variables of the monic type finite Rn chain and confirm that 
the proposed algorithm computes the generalized eigenvalues of a given matrix pencil and its convergence 
speed depends on the parameters s^' and Kf +n . 

Figure [l] shows the result with the parameters s^) = 1.01 for alH > and k„ = 1 for all n > 4, where q\P 
and e„ are the variables of the monic type Rn chain, v„ are the diagonal elements of and ic„ are the 
subdiagonal elements of (see equations | |3.6b) and | 3.6c) ). We can confirm that V„ and Ztf„ converge 
linearly to the eigenvalues and zero, respectively. Since the shift parameter s^' is not so close to the minimal 
eigenvalue 6/5 = 1.2, the stopping criterion is satisfied at t = 4605. 

Figure [2] shows the result with more suitable parameters: = 1.19 for all t > and k„ = -10000 for all 
n > 4. The convergence speed is much faster than the former example; the stopping criterion is satisfied at 
t = 48. Table[T]shows the computed eigenvalues. 

Example 2. Next, the test cases for N = 512, 1024, 2048, 4096, 8192 were computed by two methods. By 
these examples, we will compare the computation time and the accuracy of the proposed algorithm with a 
routine called DSYGV in LAPACK 3.4.2 Il28l . DSYGV computes the generalized eigenvalues of a given matrix 
pencil (A, B) in double precision, where A is real symmetric and B is real symmetric and positive definite. 
Internally, DSYGV computes the Cholesky factorization B = LL T , where L is a lower triangular matrix, trans- 
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Table 1: The eigenvalues computed by Algorithm]!] The parameters are 



1.19 and k„ = -10000 for all t > and n > 4. 



Computed eigenvalues True eigenvalues 



1.9999999999999998 
1.4999999999999991 
1.3333333333333335 
1.2500000000000000 
1.2000000000000000 



2.0000000000000000 
1.5000000000000000 
1.3333333333333333 
1.2500000000000000 
1.2000000000000000 
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Table 2: The results of the computation by Algorithm[T]for the generalized eigenvalue problems of (K^ + 2%, + 1^ ) . 



Problem size (N) 


512 


1024 


2048 


4096 


8192 


Computation time [sec] 
Maximum relative error 
Average relative error 


0.0958 
3.109 xlO" 15 
1.344 xlO" 16 


0.392 
3.405 x 10" 15 
1.211 x 10" 16 


1.58 

1.776 xHT 15 
1.154 xlO" 16 


6.24 
3.701 x 10" 15 
1.072 xlO" 16 


24.6 
2.043 x 10" 14 
1.129 xlO" 16 



Table 3: The results of the computation by DSYGV in LAPACK for the generalized eigenvalue problems of (K^ + 2I^,K^ + In). 



Problem size (N) 


512 


1024 


2048 


4096 


8192 


Computation time [sec] 
Maximum relative error 
Average relative error 


0.162 
3.664 x 10" 15 
6.673 x 10" 16 


1.92 
6.815 x 10" 15 
8.469 x 10" 16 


30.3 
1.304 x 10" 14 
1.035 xHT 15 


307 
1.684 x 10" 14 
1.276 x 10" 15 


2400 
2.949 x 10" 14 
1.508 x 10" 15 



forms the generalized eigenvalue problem A<p = xBq> to the eigenvalue problem L -1 AL~ J (L T ep) = x(L T q>) 
and solves the eigenvalue problem. We should remark that, even if A and B are both tridiagonal, L -1 AL~ T is 
a dense matrix in general. Hence, we expect that DSYGV spends much time for large problems. On the other 
hand, the proposed algorithm preserves the tridiagonal form of the matrices A^' and B^h The proposed 
algorithm will thus compute the generalized eigenvalues of tridiagonal matrix pencils fast and accurately 
for large problems. 

Tables [2] and [3] show the results of the computation by the proposed algorithm and DSYGV, respectively. 
The parameters for the proposed algorithm are s^' = (N + 2)/(N + 1) for alH > and K n = -10000 for all 
n > N - 1. In all the cases, the proposed algorithm is faster and more accurate than DSYGV. In particular, the 
proposed algorithm has an advantage in computation time for large problems. Remark that the techniques 

called deflation and splitting (if \w„ \ and |A n io„ | become sufficiently small for some n at a time t, then 
the problem can be deflated or split into two problems) were not implemented in the program used for the 
experiments. These techniques may further accelerate the proposed algorithm. 



5. Conclusion 

In this paper, we have studied the monic type Rn chain in detail and proposed a generalized eigenvalue 
algorithm for tridiagonal matrix pencils based on a subtraction-free form of the monic type finite Rn chain. 
It has been shown that, similarly to the dqds algorithm, the parameter s^' in the monic type finite Rn 
chain plays the role of the origin shifts to accelerate convergence and the proposed algorithm computes the 
generalized eigenvalues of tridiagonal matrix pencils fast and accurately. 

In ExampleEl the shift parameter is chosen ideally and all the conditions | |4.6| are satisfied. However, 
it is difficult to make this situation in general. Further improvements are thus required for practical use. First, 



in general, the condition for positivity 1 4.6 1 is not sufficient for applications; the condition does not provide 
concrete ways to choose the parameters for general cases. Second, for applying the proposed algorithm to 
general (not tridiagonal) matrix pencils, a preconditioning called simultaneous tridiagonalization (see, e.g., 
Il29l [30 j) is required. In addition to the improvements, comparisons with traditional methods should be 
discussed. These are left for future work. 
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